%% ========================================================================
%          Benefit Cost Analysis
%%=========================================================================

% Set up
run('setup.m');

% Data files
filename_plants = [ data 'structural_model_input.txt'];
filename_bids   = [ data 'bids.txt'];
filename_period = [ data 'period_statistics.txt' ];

% Read in plant and bid data
edat = etsData( filename_plants, filename_bids, filename_period );
display(edat.period);

%% Estimate ETS abatement model in treatment

% Initialize estimation object
ets = etsEstim( edat.plant, edat.bid );

% Plant random effects
ets.modnum = 1;
ets = ets.estimate;
display(ets.fits{1});

% Plant fixed effects
ets.modnum = 2;
ets = ets.estimate;
display(ets.fits{2});
 
% Plant-period fixed effects
ets.modnum = 3;
ets = ets.estimate;
display(ets.fits{3});

% Plant-period fixed effects + Hetero
ets.modnum = 4;
ets = ets.estimate;
display(ets.fits{4});

% Model 5: plain plant-period fixed effects
ets.modnum = 5;
ets = ets.estimate;
display(ets.fits{5});

ets = ets.estimateControlStringency;

%% Benefit-Cost Analysis of Scaled-up ETS in Surat
status_quo_cap = 240000;

% Run counterfactuals for table
pct_reductions = [0.1, 0.3, 0.5];

% if you want to calculate marginal costs
%pct_reductions = [-2/240 2/240]
cap = [ status_quo_cap status_quo_cap*(1-pct_reductions) ];

ctr = etsCounter( ets, cap );

ctr.scenario.regime = 'Market';
ctrM = ctr.run;  

ctr.scenario.regime  = 'Command';
ctr.scenario.ermodel = 'heat_eps';
ctrC4 = ctr.run;

counterfactuals = { ctrM, ctrC4 };

outTable = counterfactuals{1}.output_tk;
outTable.cid = ones(height(outTable),1);
% Aggregate emissions within scenario-period-emissions levels
outTableAgg = grpstats(outTable,{'cid','id_level'},{'mean'},...
    'datavars',{'Et','Zt','Pt'});
outTableAgg = sortrows(outTableAgg,{'cid','id_level'},{'ascend','ascend'});
outTableAgg.Zt = round(outTableAgg.mean_Zt / 1e6,2);
ets_costs = outTableAgg.Zt;

outTable = counterfactuals{2}.output_tk;
outTable.cid = ones(height(outTable),1);
% Aggregate emissions within scenario-period-emissions levels
outTableAgg = grpstats(outTable,{'cid','id_level'},{'mean'},...
    'datavars',{'Et','Zt','Pt'});
outTableAgg = sortrows(outTableAgg,{'cid','id_level'},{'ascend','ascend'});
outTableAgg.Zt = round(outTableAgg.mean_Zt / 1e6,2);
cc_cost = outTableAgg.Zt(end);

%% Benefit cost calculation

cost_savings = dictionary(["ets", "cyc_avg", "cyc_var", "bf_avg", "bf_var", "scr_avg", "scr_var",  "esp_avg", "esp_var"], ...
    [0, 5.60, 1.69, 8, 2.41, 10.99, 3.31, 69.09, 20.82]);

for abatement_type = keys(cost_savings)'
    abatement_type

    value_mat = zeros(20, 3); %matrix for storing results

    for col_num = 1:3

        pct_reduction = pct_reductions(col_num);
        fprintf("Emissions reduction: %d%%\n", pct_reduction*100)

        ets_cost = ets_costs(end-col_num);

        % Panel A
        fprintf("==== Panel A: Annual costs from Scale-up of Emissions Trading ====\n")
        INR2USD = 70;
        num_plant = 156;
        monitoring_costs = 5000;

        if abatement_type == 'ets'
             abatement_cost_change = (ets_cost-cc_cost) * 12 / num_plant / INR2USD * 1e6;
        else
            abatement_cost_change = (status_quo_cap*pct_reduction*cost_savings(abatement_type)) * 12 / num_plant / INR2USD;
        end
       
        total_costs = (monitoring_costs + abatement_cost_change) * 906;
        total_costs
        fprintf("1. Monitoring costs (per plant): %.0f USD/year\n", monitoring_costs)
        value_mat(1, col_num) = monitoring_costs;
        fprintf("2. Abatement cost change (per plant): %.0f USD/year\n", abatement_cost_change)
        value_mat(2, col_num) = round(abatement_cost_change, 0);
        fprintf("3. Total costs: %.2f million USD/year\n", total_costs / 1e6)
        value_mat(3, col_num) = round(total_costs / 1e6,2);

        % Panel B
        fprintf("==== Panel B: Reduction in pollution ====\n")
        ambient_pm_conc = 88.5;
        industry_share = 0.32;
        reduction_ambient_pm = ambient_pm_conc * industry_share * pct_reduction;
        fprintf("1. Ambient PM2.5 concentration: %.1fµg/m^3\n", ambient_pm_conc)
        value_mat(4, col_num) = round(ambient_pm_conc,1);
        fprintf("2. Industry share: %.2f\n", industry_share)
        value_mat(5, col_num) = round(industry_share,2);
        fprintf("3. Reduction in industry PM2.5: %.0f%%\n", pct_reduction*100)
        value_mat(6, col_num) = round(pct_reduction*100,0);
        fprintf("4. Reduction in ambient PM2.5: %.2fµg/m^3\n", reduction_ambient_pm)
        value_mat(7, col_num) = round(reduction_ambient_pm,2);

        % Panel C
        fprintf("==== Panel C: Population and value of statistical life ====\n")
        population = 7.5*1e6;
        life_expectancy = 70;
        vsl = 665000;
        value_of_1year = vsl / life_expectancy;
        fprintf("1. Population: %.1f million\n", population/1e6)
        value_mat(8, col_num) = round(population/1e6,1);
        fprintf("2. Life expectancy: %.0f years\n", life_expectancy)
        value_mat(9, col_num) = round(life_expectancy,0);
        fprintf("3. Value of statistical life: %.0f USD\n", vsl)
        value_mat(10, col_num) = round(vsl,0);
        fprintf("4. Value of 1 year of life: %.0f USD/year\n", value_of_1year)
        value_mat(11, col_num) = round(value_of_1year,0);

        % Panel D
        fprintf("==== Panel D: Gain in life-years and valuation of benefits ====\n")
        mortality_impact = 0.098;
        gain_in_life_expectancy = mortality_impact * reduction_ambient_pm;
        gain_in_life_exp_per_year_ets = gain_in_life_expectancy / life_expectancy;
        total_gain_life_years = gain_in_life_exp_per_year_ets * population;
        value_gain_life_years = total_gain_life_years * value_of_1year;
        fprintf("1. Mortality impact: %.3f years/(µg/m^3)\n", mortality_impact)
        value_mat(12, col_num) = round(mortality_impact,3);
        fprintf("2. Gain in life expectancy: %.3f years\n", gain_in_life_expectancy)
        value_mat(13, col_num) = round(gain_in_life_expectancy,3);
        fprintf("3. Per year of ETS: %.3f years\n", gain_in_life_exp_per_year_ets)
        value_mat(14, col_num) = round(gain_in_life_exp_per_year_ets,3);
        fprintf("4. Total gain in life-years: %.0f years\n", total_gain_life_years)
        value_mat(15, col_num) = round(total_gain_life_years,0);
        fprintf("5. Value of gain in life-years: %.0f million USD\n", value_gain_life_years/1e6)
        value_mat(16, col_num) = round(value_gain_life_years/1e6,0);

        % Panel E
        fprintf("==== Panel E: Benefit-cost ratio ====\n")
        benefit_cost_ratio = value_gain_life_years/total_costs;
        fprintf("1. Ebenstein et al. (2017): %.0f:1\n", benefit_cost_ratio)
        value_mat(17, col_num) = round(benefit_cost_ratio,0);
        literature = ["Pope et al. (2009)", "Correia et al. (2013) ", "Apte et al. (2018)"];
        i = 1;
        for mortality_impact = [0.061, 0.035, 0.012]
            gain_in_life_expectancy = mortality_impact * reduction_ambient_pm;
            gain_in_life_exp_per_year_ets = gain_in_life_expectancy / life_expectancy;
            total_gain_life_years = gain_in_life_exp_per_year_ets * population;
            value_gain_life_years = total_gain_life_years * value_of_1year;
            benefit_cost_ratio = value_gain_life_years/total_costs;
            fprintf("%.0f. %s: %.0f:1\n", i+1, literature(i), benefit_cost_ratio)
            value_mat(17+i, col_num) = round(benefit_cost_ratio,0);
            i = i + 1;
        end
        fprintf("\n")

    end

    %% Format strings

    jf = java.text.DecimalFormat;
    str_cell = cell(20,3);
    for i=1:20
        for j=1:3
            value = value_mat(i,j);
            if (i == 14)
                str_cell{i,j} = num2str(value,'%.3f');
            elseif abs(value) >= 1000
                str_cell{i,j} = char(jf.format(value));
            elseif (abs(value) > 1) && (abs(value) < 20)
                str_cell{i,j} = num2str(value,'%.1f');
            elseif (abs(value) < 1)
                str_cell{i,j} = num2str(value,'%.1f');
            else
                str_cell{i,j} = num2str(value);
            end
        end
    end
    disp(str_cell)


    %% Write LaTeX Table
    if abatement_type == 'ets'
        filename = [ tables '/Table_6.tex'];

        % Open table file
        fid = fopen(filename, 'w');
        % Header
        fprintf(fid, '\\begin{tabular}{lcccll}\n');
        fprintf(fid, '\\toprule\n');
        fprintf(fid, '  & \\multicolumn{3}{c}{Emissions reduction} & & \\\\\n');
        fprintf(fid, '\\cmidrule(lr){2-4}\n');
        fprintf(fid, '  & 10\\%% & 30\\%% & 50\\%% & Units & Source \\\\\n');
        fprintf(fid, '  & (1) & (2) & (3) & (4) & (5)  \\\\ \\midrule \\addlinespace\n');
    
        % Panel A
        fprintf(fid, '\\multicolumn{6}{c}{\\textit{Panel A. Annual costs from scale-up of emissions trading}} \\\\ \\addlinespace\n');
        fprintf(fid, '1. Monitoring costs per plant & $%1$s$ & $%1$s$ & $%1$s$ & \\$/year & \\multicolumn{1}{l}{Author estimates}\\\\\n', str_cell{1,1});
        fprintf(fid, '2. Abatement cost $\\Delta$ per plant & $%s$ & $%s$ & $%s$ & \\$/year & Author estimates \\\\\n', str_cell{2,:});
        fprintf(fid, '3. Total costs & $%s$ & $%s$ & $%s$ & \\$m/year & \\multicolumn{1}{l}{= (A1 + A2)$\\times$906} \\\\\n', str_cell{3,:});
        fprintf(fid, '\\addlinespace\n');
    
        % Panel B
        fprintf(fid, '\\multicolumn{6}{c}{\\textit{Panel B. Reduction in pollution}} \\\\ \\addlinespace\n');
        fprintf(fid, '1. Ambient PM$_{2.5}$ conc. & $%1$s$ & $%1$s$ & $%1$s$ & $\\mu g / m^3$ &  \\multicolumn{1}{l}{Guttikunda et al. (2019)}\\\\\n', str_cell{4,1});
        fprintf(fid, '2. Industry share & $%1$s$ & $%1$s$ & $%1$s$ & & \\multicolumn{1}{l}{Guttikunda et al. (2019)} \\\\\n', str_cell{5,1});
        fprintf(fid, '3. Reduction in industry PM$_{2.5}$ & $%s$ & $%s$ & $%s$ & \\%% & \\\\\n', str_cell{6,:});
        fprintf(fid, '4. Reduction in ambient PM$_{2.5}$ & $%s$ & $%s$ & $%s$ & $\\mu g / m^3$ & \\multicolumn{1}{l}{= B1$\\times$B2$\\times$B3} \\\\\n', str_cell{7,:});
        fprintf(fid, '\\addlinespace\n');
    
        % Panel C
        fprintf(fid, '\\multicolumn{6}{c}{\\textit{Panel C. Gain in life-years}} \\\\ \\addlinespace\n');
        fprintf(fid, '1. Mortality impact & $%1$s$ & $%1$s$ & $%1$s$ & $\\nicefrac{\\text{years}}{(\\mu g / m^3)}$ & Ebenstein et al. (2017) \\\\\n', str_cell{12,1});
        fprintf(fid, '2. Gain in life expectancy & $%s$ & $%s$ & $%s$ & years & = C1$\\times$B4 \\\\\n', str_cell{13,:});
        fprintf(fid, '3. Life expectancy & $%1$s$ & $%1$s$ & $%1$s$ & years & \\multicolumn{1}{l}{World Bank} \\\\\n', str_cell{9,1});
        fprintf(fid, '4. Per year of ETS & $%s$ & $%s$ & $%s$ & years & = C2 / C3 \\\\\n', str_cell{14,:});
        fprintf(fid, '5. Population & $%1$s$ & $%1$s$ & $%1$s$ & m & \\multicolumn{1}{l}{World Pop. Review} \\\\\n', str_cell{8,1});
        fprintf(fid, '6. Total gain in life-years & $%s$ & $%s$ & $%s$ & years & = C4$\\times$C5 \\\\\n', str_cell{15,:});
        fprintf(fid, '\\addlinespace\n');
        
        % Panel D
        fprintf(fid, '\\multicolumn{6}{c}{\\textit{Panel D. Value of gain in life-years}} \\\\ \\addlinespace\n');
        fprintf(fid, '1. Value of statistical life & $%1$s$ & $%1$s$ & $%1$s$ & \\$ & \\multicolumn{1}{l}{Nair et al. (2021)}\\\\\n', str_cell{10,1});
        fprintf(fid, '2. Value of 1 year of life & $%1$s$ & $%1$s$ & $%1$s$ & \\$/year & \\multicolumn{1}{l}{= D1 / C3}\\\\\n', str_cell{11,1});
        fprintf(fid, '3. Value of gain in life-years & $%s$ & $%s$ & $%s$ & \\$m & = C6$\\times$D2 \\\\\n', str_cell{16,:});
        fprintf(fid, '\\addlinespace\n');
    
        % Panel E
        fprintf(fid, '\\multicolumn{6}{c}{\\textit{Panel E. Benefit-cost ratio}} \\\\ \\addlinespace\n');
        fprintf(fid, '1. Ebenstein et al. (2017) & $%s$:1 & $%s$:1 & $%s$:1 & & = D5 / A3\\\\\n', str_cell{17,:});
        fprintf(fid, '2. Pope, Ezzati and Dockery (2009) & $%s$:1 & $%s$:1 & $%s$:1 \\\\\n', str_cell{18,:});
        fprintf(fid, '3. Correia et al. (2013) & $%s$:1 & $%s$:1 & $%s$:1 \\\\\n', str_cell{19,:});
        fprintf(fid, '4. Apte et al. (2018) & $%s$:1 & $%s$:1 & $%s$:1 \\\\\n', str_cell{20,:});
        fprintf(fid, '\\addlinespace\n');
    
        % Footer
        fprintf(fid, '\\bottomrule\n');
        fprintf(fid, '\\end{tabular}\n');
        fclose(fid)
    end
end
